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Abstract 

Dirac-S distributions are often crucial components of the solid-fluid coupling operators 
in immersed solution methods for fluid-structure interaction (FSI) problems. This is 
certainly so for methods like the Immersed Boundary Method (IBM) or the Immersed 
Finite Element Method (IFEM) , where Dirac-5 distributions are approximated via smooth 
functions. By contrast, a truly variational formulation of immersed methods does not 
rcqiiire the use of Dirac-(5 distributions, either formally or practically. This has been 
shown in the Finite Element Immersed Boundary Method (FEIBM), where the variational 
structure of the problem is exploited to avoid Dirac-J distributions at both the continuous 
and the discrete level. 

In this paper, we generalize the FEIBM to the case where an incompressible Newtonian 
fluid interacts with a general hyperelastic solid. Specifically, we allow (i) the mass 
density to be different in the solid and the fluid, (ii) the solid to be either viscoelastic 
of differential type or purely elastic, and {Hi) the solid to be either compressible or 
incompressible. At the continuous level, our variational formulation combines the natural 
stability estimates of the fluid and elasticity problems. In immersed methods, such 
stability estimates do not transfer to the discrete level automatically due to the non- 
matching nature of the finite dimensional spaces involved in the discretization. After 
presenting our general mathematical framework for the solution of FSI problems, we 
focus in detail on the construction of natural interpolation operators between the fluid 
and the solid discrete spaces, which guarantee semi-discrete stability estimates and strong 
consistency of our spatial discretization. 

Keywords: Fluid Strueutrc; Interaction; Immersed Boundary Methods; Immersed 
Finite Element Method; Finite Element Immersed Boundary Method 



1. Introduction 

There are several approaches to the solution of general fluid-structure interaction 
(FSI) problems. Among these we find a set of methods, which we will call immersed 
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methods, for which the discretization of the fluid domain is completely independent of 
that of the solid. These methods are more recent than and stand in contrast to established 
methods like the arbitrary Lagrangian-Eulerian (ALE) ones (see 



e-g-, 



Hughes et al. 



1981 ), where the topologies of the solution grids for the fluid and the solid are constrained. 
When the flow can be modeled as a linear Stokes flow, reduced methods, such as the 
boundary element method (see, e.g., Alouges et al. 2008 20111 can be very efficient. 



However, for more general situations immersed methods offer appealing features. 

In immersed methods the solid domain is surrounded by the fluid. When the fluid 
and solid do not slip relative to one another, these methods have three basic features: 

1. The support of the equations of motion of the fluicQ is extended to the union of 
the physical fluid and solid domains. 

2. The equations of motion of the fluid have terms that, from a continuum mechanics 
viewpoint, are body forces "informing" the fluid of its interaction with the solid. 

3. The velocity held of the immersed solid is identified with the restriction to the solid 
domain of the velocity field in the equations of motion of the fluid. 

In many respects, immersed methods can be distinguished from one another depending 
on how these three elements are treated theoretically and/or are implemented practically. 

Immersed methods were pioneered by Peskin and his co-workers (Peskin 1977 see 
also Peskin 2002 for a comprehensive account) who proposed an approach known as 
the immersed boundary method (IBM). In the IBM, the approximate solution of the 
extended fluid flow problem is obtained via a finite difference (FD) scheme. The body 
forces expressing the FSI are determined by modeling the solid body as a network of elastic 
fibers with a contractile element. As such, this system of forces has singular support (the 
boundary in the method's name) and is implemented via Dirac-i^ distributions. From a 
numerical viewpoint, the configuration of the fiber network is identified with that of a 
discrete set of points. The motion of these points is then related to the motion of the fluid 
again via Dirac-(5 distributions. Hence, Peskin's approach relies on Dirac-(5 distributions 
twice: first for the determination of the FSI force system, and again for the determination 
of the motion of the fiber network. In the method's numerical implementation, the 
Dirac-(5 distributions are aproximated as functions. Both the use of FD schemes and the 
approximation of the Dirac-^ distribution yield inconveniences that can be avoided by 
reformulating the problem in variational form and adopting corresponding approximation 
schemes such as finite element methods (FEM). 

The replacement of the FD scheme with a finite element method (FEM) was first 



proposed, almost simultaneously, bylBoffi and Gastaldi (20031, Wang and Liu (20041 



and Zhang et al. (2004). Boffi and Gastaldi (2003) show that a variational formulation of 



the problem presented by Peskin ( 1977) does not necessitate the approximation of Dirac- 
S distributions. The explicit presence of Dirac-(5 distributions pertaining to the body 
force system "disappears" naturally in the weak formulation. As far as the motion of the 
solid is concerned, the use of the Dirac-(5 distribution is unnecessary because the finite 
element solution for the fiuid velocity field can be evaluated over the solid domain on a 



pointwise basis. The thrust of the work by Wang and Liu (2004) and Zhang et al. (2004) 



^These equations are typically taken to be the Navier-Stokes equations (see, e.g., |Peskin| |1977| |Boffi| 
[and Gastaldi, 2003 Bo fB et al.[ [2008 Wang ~nd Liu[[2004[ |. However, formulations in which the fluid is 
modeled as slightly compressible have also been proposed ( |Liu et aT] |2007| |Wang et al.[ |2009| . 
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was instead that of removing the requirement that the immersed solid be a "boundary." 
The methods proposed in Wang and Liu (j2004j) and Zhang et al. (20041 apply to solid 
bodies of general topological and constitutive characteristics. However, these approaches 
still require an approximation of the Dirac-(5 distribution as a function. Specifically, 
Wang and Liu ( |2004 ) and Zhang et al. (20041 rely on the reproducing kernel particle 



method (RKPM) to approximate the Dirac-5 distribution both for the expression of the 
interaction forces and for the determination of the velocity of the immersed solid. For 



future reference, we point out that the work by Wang and Liu (2004) and Zhang et al. 



( 2004 1 pertains to systems consisting of a nearly incompressible solid body immersed in 
a Newtonian fluid. 



The generalization of the approach proposed by BofR and Gastaldi ( 2003 1 to include 
regular solid bodies, as opposed to boundaries, has been presented in various publica- 



tions culminating in the works by Heltai (2006), Bofln et al. (20081, and Heltai (2008) 



(see bibliographic references in these publications for details). The constitutive behavior 
of the immersed solid is assumed to be visco-elastic with the viscous component of the 
solid stress response being identical to that of the fluid. Another restrictive assumption 
of this work was the assumption that the fluid and the solid have the same mass den- 
sity distribution. From the viewpoint of the treatment of the interaction forces and of 
the velocity equation for the solid body, these works show that the FEM allows one to 
completely avoid dealing with Dirac-(5 distributions and their approximation. They also 
show that the velocity fleld of the solid domain can be determined variationally, i.e., in a 
way that is consistent with the FEM as a whole. However, the strength of this idea is not 
fully demonstrated by Boffi et al. (2008). This is because they choose a flnite element 
discretization of the solid domain for which the motion of the solid is determined by a 
direct evaluation the fluid's velocity at the vertices of the grid supporting the discretiza- 
tion in question. The advantage of a fully variational formulation for immersed methods 
is the fact the FEM machinery offering transparent stability results and error estimates 
becomes readily available along with the machinery developed for adaptivity. 

A fully variational formulation of the immersed problem that docs not rely on any 
approximation of the Dirac-(5 distribution has also been formally discussed by |Liu et al.| 
( [2007| ) (see E q. (40) on p. 215). However, it is not clear whether or not the variational 



formalism of Liu et al. (2007) has been implemented in actual calculations. Another 



fully variational formulation of an immersed method has been proposed by [Blanco et al.| 
(2008a). This formulation can also cope with a variety of constitutive assumptions for 
the fluid and solid domains. While some numerical results have been published for the 



case of solid structures having incompatible kinematic assumptions (see Blanco et al. 



2008b), no numerical results seems to have been published for FSI problems. 



Here we present the generalization of the approach discussed by |Bofii et al. (2008) so 
as to be applicable to the case of general visco-elastic compressible and incompressible 
bodies immersed in an incompressible fluid. The proposed scheme produces a discretiza- 
tion which is strongly consistent and stable, and can easily be extended to the case in 
which the fluid is also compressible. 

As mentioned earlier, in immersed methods the velocity of the solid is set equal to 
the restriction to the solid domain of the velocity of the fluid. When enforced variation- 
ally, this equality is weakened and transport theorems underlying the classical energy 
estimates typically obtained in continuum mechanics do not hold any longer. While 
the classical transport theorems cannot be invoked directly, we show that energy esti- 
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mates and corresponding stability results can be obtained for our proposed abstract weak 
formulation that are formally identical to the classical ones from continuum mechanics. 

The proposed variational formulation produces a natural discretization scheme which 
differs from the one presented by Boffi et al. ( j2008J in the determination of the motion 
of the solid. In BofR et al. (2008) the velocity field is evaluated at the discrete level 
on the vertices of the solid mesh. This procedure renders semi-discrete and discrete 
stability estimates nontrivial for general approximating spaces of the solid displacement. 
By contrast, the stability results we prove in the abstract weak formulation are inherited 
naturally by the discretization scheme, provided that conforming approximating spaces 
are used for the velocity and displacement fields, thus removing the assumptions on the 
the triangulation of the solid that were present in BofB et al. ( 2008 ) . 

Another original contribution of our formulation is that the treatment of the case 
of a compressible solid in an incompressible fluid is not taken as a limit case of some 
other set of constitutive assumptions. This is an important detail in that, again to the 
best of the authors' knowledge, other approaches have dealt with solid/fluid kinematical 
incompatibilities indirectly, i.e., as limit cases corresponding to some particular value of 
a tunable parameter. 

In Section [2j we present a formulation of the equations of motion of an immersed 
solid in the context of classical continuum mechanics (i.e., under strong smoothness 
assumptions) . We also offer a concise exposition of the transport theorems and associated 
energy estimates that are valid in the aforementioned classical context. In Section [3] we 
reformulate the problem in variational form and present a discussion of the formulation's 
underlying functional setting. We then prove the that proposed formulation is stable. 
In Section [4] we present the discrete formulation we derive from the proposed abstract 
variational formulation and show that the discrete formulation is strongly consistent and 
inherits the stability of the abstract formulation. Some numerical results are presented 
in Section [5] 



2. Classical Formulation 

2.1. Basic notation and governing equations 

Referring to Fig. [l] Bt represents the configuration of a solid body S§ at time t. As a 














y p 




^^p(t) 



Figure 1: Current configuration Bt of a body SS immersed in a fluid occupying the domain f2. 
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point set, Bf is a (possibly multiply connected) proper subset of a fixed control volume 
f2. The domain fl \ Bt is occupied by a fluid and we refer to Bt as the immersed body. 
The boundaries of and Bt, with outer unit normals m and n, respectively, will be 
denoted by d^l and OBt- For convenience, we select a configuration of ^ as a reference 
configuration and we denote it by -B. Both B and Bt are viewed as submanifolds of a same 
Euclidian manifold S"^, of dimension d equal to 2 or 3, covered by a single rectangular 
Cartesian coordinate system with origin at O. We denote the position of points of in 
B by s, whereas we denote the position at time i of a generic point P G ft hy xpit). 
A motion of ^ is a diffeomorphism C : B ^ Bt, x = C{s,t), with s E B, x E il, and 
t G [0,r), with T a positive real number. 

We denote by p{x,t) the spatial (or Eulerian) description of the mass density at the 
location x at time t. The function p can be discontinuous across dBt- The local form of 
the balance of mass requires that, G (0,T), 

p + pdivM = 0, X En\{dn[JdBt), (1) 

where u{x,t) = 9C,{s ,t) / dt \ is the spatial description of the material velocity 

field, a dot over a quantity denotes the material time derivative of that quantityj^ and 
where 'div' represents the divergence operator with respect to x. 

We denote by T{x,t) the spatial description of the Cauchy stress. The local form of 
the momentum balance laws require that, Vi € (0, T), T = T""" (the superscript T denotes 
the transpose) and 

divT + pb = ptt, X EVL\{dnyjdBt), (2) 

where b{x,t) describes the external force density per unit mass acting on the system. 

In addition to Eqs. ([l]) and ([2]), we also require the satisfaction of some continuity 
conditions across dBt- Specifically, we demand that the velocity field be continuous 
(corresponding to a no slip condition between solid and fluid) and that the jump condition 
of the balance of linear momentum be satisfied across dBt- For all t £ (0, T), these two 
conditions can be expressed as follows: 

uix'^ ,t) = u{x^ ,t) and 'J{x'^ ,t)n = l'{x^ ,t)n, x E dBt, (3) 

where the superscripts — and + denote limits as a; — > ac from within and without Bt, 
respectively. 

We denote by dVlo and 917 jv the subsets of dVl where Dirichlet and Neumann bound- 
ary data are prescribed, respectively. The domains OVLd and 80. are such that 

dVL = dnoy^dVLN and 9171)0 9^2^ = 0. (4) 



^In continuum mechanics, a physical body is assumed to consist of material points. Each of these 
is considered an individual entity whose position is a function of time and at which physical quantities, 
such as mass density or momentum density, can be defined. By definition, the material time derivative 
of a property of a material point (e.g., momentum), is the time rate of change of that property measured 
while following the material point in question. The material time derivative of a (scalar-, vector-, or 
tensor-valued) field of the type 4' = <^{s,i), with s S B, is simply </> = d4>/dt. In the case of a scalar- 
valued function ij} = ^{x,t), with x Q, ip = dip/dt + (grad^/)) ■ u, where 'grad' is the gradient with 
respect to x, u{x,t) is the (material) velocity field, and '■' denotes the standard inner product for vectors 
fields. For a vector-valued function w{x,t), w = dw/dt + (gradK;)tt, where '(gradio)tt' denotes the 
action of the second order tensor gradio on the velocity field it. 
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We denote by Ug{x,t), with x G d^jj, and by Tg{x,t), with x € dftN, the prescribed 
values of velocity (Dirichlet data) and traction (Neumann data), respectively, i.e.. 



u{x,t) — Ug{x,t), ioT X G dflD, and T{x,t)m{x,t) — Tg{x,t), for x E OUn, 

(5) 

where the subscript g stands for 'given.' 

Using the principle of virtual work and letting v denote any admissible variation of 
the field u, Eqs. ([2| and ^ can be written as follows: 



/ p{u — b) ■ vdv + / T • grad vdv — Tg ■ v da = 0, 

Jn Jn JdnN 



(6) 



where da and dv represent infinitesimal area and volume elements, respectively. We can 
reformulate Eq. ([l]) in variational form as follows: 



P_ 
n\P 



divu \qdv — 0, 



(7) 



where, from a physical viewpoint, q represents any admissible variation of the pressure 
in the system. In the case of incompressible materials, p = and Eq. ([7| yields the 
traditional weak form of the incompressibility constraint, namely, /^^ qdivudv — 0. 

2.2. Constitutive behavior 

2.2.1. Constitutive response of the fluid. 

We assume that the fluid is linear viscous and incompressible with uniform mass 
density pf. Denoting by Tf the constitutive response function of the Cauchy stress of the 
fluid, we have (see, e.g., Gurtin et al., 2010) 

tf = -pl + 2MtD, D-i(L + LT), (8) 

where p is the pressure of the fluid, I is the identity tensor, /if > is the dynamic 
viscosity of the fluid, and L = gradw, and where a "hat" (T) is used to distinguish the 
constitutive response function for T from T itself. For convenience, we denote by Tj the 
viscous component of Tf, i.e., 

t- = 2/if D = /if (L+LT). (9) 
Incompressibility demands that pf = so that Eq. ([T]) yields the kinematic constraint 

divM = for a; e 17 \ St. (10) 



Under these conditions, p is a Lagrange multiplier allowing us to enforce Eq. (10). In 
addition, Eq. (10) also implies that trL — so that the term Tf in Eqs. ([8) is the 
deviatoric part of the Cauchy stress in the fluid. 
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2. 2. 2. Constitutive response of the solid. 

We assume that the body is viscoelastic of differential type. The response function 
for the Cauchy stress of the sohd is assumed to have the following form: 



(11) 



where T| and denote the elastic and viscous parts of Tg, respectively. The viscous 
part of the behavior is assumed to be of the same type as that of the fluid, that is, 



^2^l,D = /i, (L + 



(12) 



where /Zg > is the dynamic viscosity of the solid. We do include the possibility that 
/is might be equal to zero, in which case the solid behaves in a purely elastic manner. 
As far as T| is concerned, we assume that it is given by a strain energy potential. To 
describe this part of the behavior in precise terms, we introduce the first Piola-KirchhofF 



stress tensor, denoted by P and defined as (see, e.g., Gurtin et al. 20101: 

P = JTF"^, 

where J = det F, and the tensor F, called the deformation gradient, is defined as 

^^ dC{s,t) 
ds 



(13) 



(14) 



As is standard in continuum mechanics (see, e.g., Gurtin et al. 2010), we require J to 
satisfy the following assumption: 



J(s, t) > J,„ > \/s e B and\/t e [0, t). 



(15) 



Therefore, F always admits an inverse, as required for Eq. ( 13 ) to be meaningful. Hence 



letting Pg — JT^F~"^ denote the constitutive response function for the elastic part of the 
first Piola-Kirchhoff stress tensor, as is typical in elasticity, we assume that there exists 
a function V[^s'^(F) such that 



df 



(16) 



where VKf is the constitutive response function of the volume density of the elastic strain 
energy of the solid. To satisfy invariance under changes of observer, must be a 
function of an objective strain measure such as C = F"'"F. In addition, if the solid is 
isotropic, will be taken to be a function of the principal invariants of C. Finally, if 
the solid is incompressible, then its stress response is determined by deformation only 
up to a hydrostatic component. In this case, the constitutive response function for the 
solid has the form 

tg = -pl+tf +tg", (17) 

where p is the Lagrange multiplier needed to enforce incompressibility, Tg is given by 



Eq. (121, and T| is obtained from Eq. (16) via Eq. (13 1 
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2.2.3. Elastic strain energy and dissipation 

While more general cases can be considered, we assume that Wjf (F) is a convex 
function over the set of second order tensor with positive determinant. As far as the 
viscous part of the behavior is concerned, we have already assumed that /if > and 
/is > 0. These conditions imply that 

L > 0, L > (18) 



for all L 7^ 0. Equations ( 18 ) imply that the viscous part of the behavior is dissipative. 



2.2.4. Mass density distribution 

As a last aspect of the formulation related to constitutive behavior, we will denote 

by 

Pso=Pso(s), seB, (19) 

the referential (or Lagrangian) description of the mass density of the solid. While Eq. ^ 
holds for the solid as well as the fluid, the local form of the balance of mass for a solid 
is typically expressed in Lagrangian form as follows: 

Pso(s) = Ps(a5,i)|^^^(^ t)J(s,t), seB, (20) 

where ^3(3^, t) is the spatial description of the mass density of the solid. We will indicate 
the general mass density of the system with p = p(x, t), with the underlying assumption 
that 

r ,N fpf' foTxen\Bt, 
p{x,t) ^ i (21) 
\^Ps[x,t), tor X e Bt, 

where, as stated earlier, pf is a constant. 

2.3. Transport theorems 

Transport theorems are kinematic results pertaining to the time differentiation of 
integrals over time-dependent domains. These results are useful in the discussion of 
energy estimates. 

Theorem 1 (Transport theorem for generic time dependent domains). Let n{t) G , 
with d = 2, 3 and t S K, be a regular, possibly multiply- connected time- dependent domain 
with boundary dQ{t) . Let m be the unit normal field orienting dfl{t), outward with respect 
to fl(t). Let V be the velocity of dVL{t) according to some convenient time-parametrization 
of dfl{t). Let (j){x,t), with x € n{t) be a .smooth field defined over n{t). Then we have 

^ '^'x,t)dv^ [ ^'^(^'^) I ,p{x,t)u -mda. (22) 



./o(t) ' Jn(t) dt Jdh(t) 



Theor em [T] is is a well-known result whose proof is available in various textbooks 



(see, e.g., Truesdell and Toupin 1960 Gurtin et al., 2010). The following form of the 



transport theorem is a simple but new result that is particularly suited for the analysis 
of the motion of immersed bodies. The proof of the theorem below can be found in the 
Appendix. 
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Theorem 2 (Transport theorem for a control volume containing an immersed domain). 
Let and Bt be the domains defined in Section \2.1\ That is, let Bt be the current 
configuration of a body immersed in the control volume il. Let ipix, t) denote the Eulerian 
description of a density per unit mass defined over D, Z) Bt, smooth over the interiors of 
fl\Bt and Bt but not necessarily continuous across dBt- Also let p{x,t) be the Eulerian 
description of the mass density distribution, which need not be continuous across dBt . 
Then 

— / p-0dw+ / pipu-mda^ i pi[)dv. (23) 
^ In Jan Jn 



dt 



Remark 1 (Generality of Theorem [2| . Theorem [2] is a straightforward but nontrivial 
result implied by the combined application of Theorem [T] and the balance of mass. One 
crucial element of Theorem [2] is that no special assumption on the behavior of the mass 
density was necessary. That is, Theorem [2] is valid whether or not p is constant or the 
fluid flow is steady. 

2.4- Theorem of power expended 

In this paper we propose an immersed method for the numerical solution of the 
problem governed by Eqs. ^ and ([7| under standard physical assumptions concerning 
the constitutive behavior of the fluid and of the immersed solid. We will discuss energy 
estimates and associated stability properties for the proposed method. To facilitate this 
discussion, it is useful to relate the power supplied to the system and the system's time 
rate of change of kinetic energy. Such a relationship is typically referred to as the theorem 



of power expended (see, e.g., Gurtin et al. 2010). Here we derive a form of the theorem 



of power expended that fits our purposes. Before doing so we introduce the following 
definitions: 



K{x,t):^lp{x,t)u'{x,t) and J-{x,t) = {J' ^ * 

I , tor a; e , 



(24) 



where k is the kinetic energy density per unit volume and v? .= u ■ u. 

Theorem 3 (Theorem of power expended for a control volume with an immersed do- 



main). Let and Bt be the domains defined in Section 2.L That is, let Bt be the current 
configuration of a body immersed in the control volume 51. Let the motion of the system 
be governed by Eqs. ^ and Q. Then 



d r , r , d 



pb-udv+ Tg-uda=--- Kdv+ KU-mda+--- dV + / T^-Ldw. 

n Jdniv Jq Jqq dt Jb Jn 

(25) 

Proof of Theorem Replacing v with u in Eq. ^ and rearranging, we obtain 

/ pb-udv+ / Tg-uda= / pu-udv+ / T-Ldw, (26) 
Jn JnN Jn Jn 

where we have used the fact that L = gradit. We now observe that 



pit ■ u — ^pu^, 
9 



(27) 



where the hne over simply denotes the fact that the material time derivative (denoted 
by the dot over the line) must be applied to the quantity under the line, namely, u^. 
Therefore, recalling that, by the first of Eqs. (24), k = \pv? , we have that 



pit ■ u dv 



pu ■ u dv 



d 
dt 



ndv 



KU ■ m da, 



(28) 



where, to obtain this last expression, we have used Theorem [2] Next, using the consti- 
tutive equations in Section 2.2 for x Cz Q \ Bt, i.e., in the fluid, we have that 



T • L = -pi • L + • L 



T - L = -pdivM + tj' • L 



T • L = t^ L, (29) 



where we have used the fact that, in the fluid, divtt = due to incompressibility. For 
X € Bt, i.e., in the solid, we would normally have to distinguish between the compressible 
and incompressible cases. However, the final result is the same due to the fact that, in 
the incompressible case, the Lagrange multiplier p does not contribute to the stress power 
as was shown in Eqs. (29). Therefore in the solid we have 

T- L 



TM + Tf • L. 



(30) 



Using Eqs. (29 1 and (30) along with the definition in the second of Eqs. (24), whether 
the solid is compressible or not, we have 



T-Ldv = 



n\Bt 



t^Ldw+ / f^-ldv 



Bt 



T-Ldv = / • Ldv 



Tt ■ Ldv 



Bt 



(31) 



• Ldv. 



Bt 



Next, recalling that L — FF ^, we recall that 



Bt 



TMdt;= / JT"-FF-My 



TMdi;= / Pf-FdK 



Bt 



(32) 



where we have used Eq. (13) along with the tensor identity A • BC = AC • B. Using 



Eq. (16) we see that • F = Wjf, so that, combining the results in the last of Eqs. (31) 
and ( 32 ) , we can write 



T-Ldw= / T"-Ldw+ / W:dV. 



(33) 



We recall that VKf = Wi{s,t), s e B, so that = dW^/dt. Therefore, observing 
that S is a fixed domain (with fixed boundary), using Theorem [l] with the identification 
Ct —i' B, we can express Eq. (33) as follows: 



T- Ldw 



T" • Ldt; + — 
dt 



W! dV. 



(34) 



The proof can be now concluded by substituting the last of Eqs. (28) and (34|) into 
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Eq. (26), which yields Eq. (25). 



□ 



Lemma 1 (Dissipation inequality) . Referring to Theorem if the system is provided 
no power input, i.e., if 

Ug =0, Tg = 0. and 6 = 0, (35) 
then, for all admissible motions of the system, 

^ [ Kdv+ [ KU- mda + [ W!^ dv < 0. (36) 
at Jn Jaa« ^t Jg 

Proof of lemma [7} Inequality ( 36 1 is a direct consequence of Theorem [3] and Eqs. 

□ 

Remark 2 (Energy estimates). Lemma [l| plays an important role in that it provides the 
form of the energy estimates and corresponding stability condition we strive to satisfy in 
the proposed numerical scheme. 



3. Abstract Variational Formulation 



We now reformulate the governing equations as a problem to be solved via a general- 
ization of the approach proposed by ^Boffi et al. (2008). For simplicity, we first consider 
the case with ^ incompressible and then the case with ^ compressible. In either case, the 
principal unknown describing the motion of the solid is the displacement field, denoted 
by w and defined as 

w{s,t) -.^ Cis,t) - s, seB. (37) 
The displacement gradient relative to the position in B is denoted by H: 



H 



dw 
ds 



H = F - I. 



Equation (37) implies 



w{s, t) — u{x, t)\ 



(38) 



(39) 



Remark 3 (Eulerian-Lagrangian information exchange and numerical approximation). 
On the one hand, u{x,t) and w{s,t) can be said to carry the same information in that 
both represent velocity. On the other, the information carried by u{x,t) and w{s,t) 
is "packaged" in fundamentally different ways in that u(x,t) is Eulerian and w{s,t) 
is Lagrangian. Equation (39) "regulates" how the information exchange occurs. As 



long as pointwise values of u{x, t) are available, given an s g i? and having a full- 
field representation of ^(s,t), then the evaluation of Eq. ([39]) is straightforward. The 
evaluation of w[s,t) is not straightforward when the field u{x,t) is not available as a 



field. As stated by Peskin (see the beginning of Section 6 in Peskin 



1977^ 



The Lagrangian mesh upon which the boundary forces /j, , and the boundary 
configuration a;^ are stored as points which do not coincide with fluid mesh 
points. We therefore have the problem of interpolating the velocity field from 
the fluid mesh to the boundary points and spreading the boundary forces 
from the boundary points to the nearly mesh points of the fluid. 



^In the cited passage, ccj, is a discrete set of points on the immersed boundary at which the forces 
/fc responsible for expressing the FSI are defined. 
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To understand the quote, it is important to recall that Peskin is solving the problem 
via FD. Therefore u{x,t) is available only as a set of discrete values at the mesh point 
defining the FD solution domain. To interpolate the discrete values of it at a point Xk 
on the immersed boundary (not coinciding with the mesh points for the FD grid), Peskin 



Peskin 19771: 



presents what is Eq. (39 1 in this paper in terms of the Dirac-(5 distribution (see Eq. (2.9) 
in 



dxk/dt — u{xk,t) = / u{x,t)5{x — Xk) dv, 



(40) 



where dx^/dt corresponds to what we would denote by ■u;(a;fe,i), and where the integral 
defines the action of the Dirac-(5 distribution on the function u{x^t). In Section 6 of 
Peskin (1977), the 5 in Eq. (40) is replaced by an actual function whose purpose is to 
approximate the behavior of the 5 and allow one to carry out the convolution integral 
explicitly. This strategy allows one to interpolate the discrete velocity field information 
at points that are not on the FD grid. What is important to notice here is that, formally. 



Eq. ( 40 ) is Eq. ( 39 1 , i.e. , they serve the same purpose of transferring Eulerian information 



into Lagrangian information. Peskin's rationale for choosing to work with Eq. ( 40 1 vs 



Eq. ( 39 ) is due to the nature of his numerical scheme. We therefore maintain that any 
numerical approximation scheme for which the fields u{x,t) and C,{x,t) are known as 
fields^ does not need to confront the issue of introducing and, a fortiori^ approximating 
Dirac-(S distributions. In this paper, the immersed problem is solved by FEM and there- 
fore it does not require the introduction of the Dirac-(5 distribution either at a formal or 
at a practical level. In our proposed approach we enforce Eq. (39) weakly, consistently 
with the variational nature of the solution method we adopt. 

3.1. Functional setting 

The principal unknowns of our fluid-structure interaction problem are the fields 

u{x,t), p{x,t), and w{s,t), with x e ft, s e B, and t e [0,T). (41) 

The functional spaces for these fields are selected as follows: 

uer = H},{nf := {u e L^inf \ V^u e L\nf'"',u\9no = ^^^^ 

pe^:^L\n), (43) 
w e ^ := e L^{Bf I S/sW e L°°(B)''^'*}, (44) 

where Va, and Vs denote the gradient operators relative to x and s, respectively. 

For convenience, we will use a prime to denote partial differentiation with respect to 
time: 

«'(.,0:.^ and «,'(.,t):.^. (45) 
Hence, in view of the discussion in the footnote on page[5j we have 

u{x,t) — u' {x,t) + (Wxu{x,t))u{x,t) and w(s,t) ^ w' {s,t). (46) 
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Remark 4 (Domains of definition of the fluid's behavior). As in every immersed method, 
a crucial element of our formulation is the extension of the domain of definition of the 



fluid's behavior to 51 as a whole. The definitions in Eqs. (42) and (43) imply that the 



fields u and p are defined everywhere in fl. Because u is defined everywhere in ft, the 
function Tj is defined everywhere in as well. For consistency, we must also extend the 
domain of definition of the mass density of the fiuid. Hence, we formally assume that 



(47) 



Remark 5 (Space of test functions for the velocity) . Referring to Eq. ( 42 1 , we will denote 



the function space containing the test functions for the velocity field by and define it 
as: 

\dxd 



H^inf := [v e L^inf I \/^v e L^inf'"', v\ann - o}. 



(48) 



Remark 6 (Functional spaces for time derivatives) . The functions u' and w' are generally 
not expected to be elements of Y and '3^, respectively. Referring to the first term in 
Eq. (|6| and the first of Eq. (46), the regularity of u' is related to the regularity of the 
given field b and of the boundary conditions. The field b is often assumed to be an 
element of H^^{^). The latter can therefore be viewed as a baseline in terms of the 
minimum regularity that u' could have. However, since the regularity of b is not the 
only factor at play, here we limit ourselves to state that u' is an element of a pivot space 
J^v such that 

f c .y/fv ^ J^v Q y* , (49) 

where J^y and are the dual spaces of and respectively. We can be more 
specific in the case of w'. We start with saying that w' is an element of a pivot space 
Ji^Y such that 

^ C C C (50) 



where and are the dual spaces of and respectively. Then, if Eq. ( 15 ) 



is satisfied, using Eq. (39) and standard Sobolev inequalities (see, e.g., Evans 2010), we 
have that, for w E 3^ and u E Y, 



^ C C H\Bf. 
In fact the H^{B)'^ norm of the displacement velocity can be controlled by 



(51) 



J B Jb 



{uo(:fdV+ / (((V,m)oC)F) dV 



= [ (ufiJoC^) ^dv+ [ (v^u{FoC')Y{JoC') ^di 

JBt JBt ^ ' 



(52) 



I + V.ti;| 



L°'{B) 



dxd VajM 



LHBt)" 



<J,„2(1 + \\\ +y sW\\l^^gy.,)\\u\\l,^s^y 

<J-,^l + \\w\\l)\\u\\'^, 



and we therefore take the pivot space Jfy to be H^{B)'^. 
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3.2. Governing equations: incompressible solid 

When the sohd is incompressible, the mass density of both the fluid and the sohd are 
constant so that p = (almost) everywhere in Q,. Cognizant of Remark |4j referring to 
Eqs. ([5]), Eqs. (42|-(|44]), and the constitutive response functions of both the fluid and 
the solid, Eqs. (6|^and ([t]) can be written as 



/ p{{u — b)-vAv+ i {ps — p{){u — b) ■ V dv 
Jn J Bf 



ff\/^vdv+ (ts - tf)-Va;t'du - / Tg-vda^O \fv e % (53) 



Bt 



and 



qdiYudv — V(7 G =S. 



(54) 



In addition to the momentum and mass balance laws, we need to enforce Eq. (39). We 
do so weakly as follows: 



w{s, t) — u{x, t) I 



y{s)dV = Q yyeJ^Y, 



(55) 



where d^ denotes the volume of an infinitesimal element of i3, and where <I>b is a constant 
with dimensions of mass over time divided by length cubed, i.e., dimensions such that, 
in 3D, the volume integral of the quantity ^gw has the same dimensions as a force. 



Remark 7 (Equation ( 55 1 and comparison with other formulations) . As discussed in the 
introduction, a key element of any fully variational formulation of immersed methods 
is (the variational formulation of) the equation enabling the tracking of the motion of 



the solid. Equation (55) is the equation in question. When discretized, it yields a set of 



ordinary differential equations (ODE) relating the degrees of freedom of the extended fluid 
domain with the degrees of freedom of the immersed domain. In practical applications, 
this relation is as general as the choice of the finite-dimensional functional subspaces 



approximating "f^ and Equation ( 55 ) plays a crucial role in ensuring that the proposed 



finite element formulation is stable. Equations similar to Eq. (55 1 have appeared in other 
variational formulation of immersed methods. With this in mind, it is important to 
remark that the set of ODE for tracking the motion of the immersed solid in 'Boffi and] 
Gastaldi (2003) was not obtained via a variational formulation. Rather, it was obtained 



by setting the value of w equal to that of u at the vertices of the triangulation discretizing 
the solid domain. The first fully variational formulation of the equations of motion of the 



solid domain was presented almost simultaneously byjHeltaij (j2006j) and |Liu et al. (2007). 
However, the present authors could not find in the literature evidence pertaining to the 



practical implementation of the the work by Liu et al. (2007). In Boffi et al. (2008) the 



solid and the fluid mass densities are the same and are equal to one (a restriction which 
was removed recently in Boffi et al. 2011). In addition, 3^ is chosen as the space of 



globally continuous piecewise affine functions over triangles in two-dimensions and over 

Under 



tetrahedrons in three dimensions (see Eq.(52) on p. 2218 in Bofln et al 



2008) 



these assumptions, Eq. (55) yields a system of ODE of the type w^it) = u{x]^,t), where 
k ranges over the index set of the vertices of the triangulation of the solid domain. That 
is, for the purpose of tracking the motion of the solid, the formulation by |Boffi et al. 
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( 2008 1 yields the same equations as those in Boffi and Gastaldi ( 2003 1 . Finally, again as 



indicated in the introduction, to the best of the authors' knowledge, the approach to the 



determination of the motion of the solid expressed via Eq. ( 55 1 has only been explicitly 
discussed by Heltai (2006), Liu et al. (2007), and Blanco et al. (2008a). However, no 



(55) 



general numerical implementations have been demonstrated. This particular aspect of 
the current formulation is one of the thrusts of this paper. A discussion of how Eq 
is practically implemented in a FEM code is presented later. 

Going back to the discussion of the problem's governing equations, we now anticipate 



that our proposed numerical approximation of Eqs. (53)-(55) is based on the use of two 



independent triangulations, namely, one of and one of B. The fields u and p, as well as 
their corresponding test functions, will be expressed via finite element spaces supported 
by the triangulation of f2. By contrast, the field w will be expressed via a finite element 
space supported by the triangulation of B. Motivated by this fact, we now reformulate 
every integral over Bt as a corresponding integral over B. Such a reformulation affects 
only Eq. ( 53 1 , which can be rewritten as 



/ p[{u — b)-vdv— / pdivvdv+ / ff-V^vdv— / Tg-vda 
Jn Jn Jn JdUN 

+ I { [Pso (s) - PiJ{^, t)] t) - b{x, t)] ■ dV 

PlF'^{s,t)-VMx)\^=^^^^t^dV ^0 Vt^e^o. (56) 



The last three terms in Eq. ( 56 1 have been written so as to explicitly express their eval- 
uation process. While it is true that, for an incompressible solid J(s, t) = 1 for sUd s £ B 
and for all t e [0, T), this occurrence may not be satisfied in an approximate formulation 
of the problem. Therefore, we prefer to retain the term J{s,t) in our formulation to 
contribute to its stability. 

Remark 8 (Dirac-i5s are not intrinsic to immersed methods). As eloquently stated by 
Boffi and Gastaldi (20031 in their introduction, "The IB method is at the same time a 



mathematical formulation and a numerical scheme." As such, and as argued in Remark[3j 
the use of the Dirac-(5 distribution was justified by convenience and a preference for a 
specific solution method rather by a necessity intrinsic to the physics of the problem. 



One of the main thrusts of the works by BofR and Gastaldi (20031; Heltai (2006); BofH 



et al. (2008); Liu et al. (2007); Blanco et al. (2008a) is precisely that of showing that 



an immersed method can be formulated without any reference whatsoever to the use of 
Dirac-(5 distributions. Again, we wish to point out that one of the objectives of the present 
work is precisely that of demonstrating an implementation technique that does not rely 
on the approximation of the Dirac-(5 distribution. This fact is one of the distinguishing 
features of our work when compared to other approaches currently in the literature (see, 
e.g.,|Wang et al.[|2009[). 



We now define various operators that will be used to state our finite element formu- 
lation. These definitions rely on the concept of duality. To make explicit the declaration 
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of the spaces in duality, we will use the following notation: 



inu] 




x,t) + 






= Ms [V^tt( 


x,t) + 




F[w] 


= \ + Vsw{s,t), 




J[w] 


= det f[w], 








_ dWiif) 








dF 


F=FH 





v'{^^^)v^ (57) 

in which, given a vector space V and its dual V* , tp and (f) are elements of the vector spaces 
V* and V, respectively, and where identifies the duality product between V* 

and V. Also, to be explicit on how certain terms depend on the selected unknown fields, 
we introduce the following shorthand notation 

(58) 
(59) 
(60) 
(61) 

(62) 

Finally, to help identify the domain and range of these operators, we establish the fol- 
lowing convention. We will use the numbers 1, 2, and 3 to identify the spaces ^, 
and respectively. We will use the Greek letter a, (3, and 7 to identify the spaces , 
and respectively. Then, a Greek letter followed by a number will identify an 
operator whose domain is the space corresponding to the number, and whose co-domain 
is in the space corresponding to the Greek letter. For example, the notations 

£a2 and £a2P (63) 

will identify a map {£0,2) from i? into Y* and the action of this map {£a2P € ^*) on the 
field p G =S, respectively. If an operators has only one subscript, that subscript identifies 
the space containing the range of the operator. For simplicity, the pivot spaces J^tfy and 
Jify and their duals will inherit the same notation as y and W. With this in mind, let 

Mai ■■ ^ , ^.{MociU,v} := piU-vdv e Jifv,'^v e %, 

Jq 

(64) 

J\fai{u) : r r* , ,{j\fai{u)w,v) := p{(ya,w)u-vdv yu^w e r.yv efo, 

Jn 

(65) 

Jq 

(66) 

Bpi-.y^ ,{BpiU,q) :=- [ qdivudv Vg e ^,Vu e r, 

Jn 

(67) 

^^r*, ,{Bjiq,u) :=- [ qdivudv Vge^,Vwer. 

Jn 

(68) 



^131 
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The operators defined in Eqs. (64|-(68) concern terms that are typical of the Navier- 
Stokes equations and will be referred to as the Navier-Stokes component of the problem. 
As in other immersed methods, these operators have their support in i7 as a whole. 

We now define those operators in our formulation that have their support over B but 
do not contain prescribed body forces or boundary terms. 



,{6Mai{w)u,v)^^ / {{ps„{s) - pfj[w])u{x) ■ v{x)}^^ 



x — s-\-w{s) 



(69) 



{5Ual{w,l,z)u,v)., 



{[{pso{s)V Mx)f-{S) 

- P{J[w\V xu{x)z{x)\ ■v{x)) 



x—s-\-'w{s) 



dV. 



x—s-\-'w{s) 



J B 



(70) 
(71) 

(72) 



We now define operators with support in B that express the coupling of the velocity 
fields defined over $7 and over B. Specifically, we have 



^.{M^3'w,y)^^:=<^B w-y{s)dV, 

J B 

M^i(w) : r ^y*, Vm G r,'iw G ^,V?/ G ^y, 



u{x) 



35 — S + TU 



y{s)dV 



Ml^{w) : .m^ r\ Vm G r,Vw G ^,Vy G ,3^Y 



.^,{M'^^,{w)y,u).^ := $b / u{x)\^^^^^^^y y{s) dV 

J B 



(73) 
(74) 
(75) 



Finally, we define the operators that express the action of prescribed body and surface 
forces. 



j"a G y\ vb G H-^{n).yTg g H-^(dnN)yv g % 



Tg ■ V da 



Pfb ■ vdv 

G T\ \fw e?V,\fbe H-\n)yv g % 
.^,{gaiw),v),^ / {ps„is)- pfj[w])b-v{x)\^^^^^^^^dz 

J B 



(76) 



(77) 
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Remark 9 (Dependence on the motion of the solid) . In defining the operators in Eqs. ( 64 1 



(77), we have used a notation meant to point out explicitly the role played by the field 



w in the evaluation of integrals over B. For the operator Ac in Eq. (72), the motion of 



the solid plays a double role, one pertaining to the elastic response of the solid (through 
w) and the other pertaining to the map (through h) functioning as a change of variables 
of integration. 

It is convenient to explicitly separate the double role of the displacement w in the 
elastic operator Aa{w,w), by reformulating it in terms of a change of variable operator 
and in terms of a purely Lagrangian elastic operator: 



Ay{w) e m;, yw e ^,Vy e 



,{A^{w),y) 



JB 



(78) 



(79) 



The operator Sa-y^h) is the map that allows us to express the duality over in terms 
of that over Y through the deformation h. As such, Sa-yih) puts into communication 
the Lagrangian and Eulerian descriptions of the motion of the immersed domain. The 



operator in Eq. ( 79 ) is a typical component of classical FEM approaches to elasticity and 
is the (fully Lagrangian form of the) stiffness operator of the immersed solid. 

One of the crucial components of any solutions method for FSI problems is the com- 
munication between the Lagrangian and Eulerian descriptions of the physics of the solid 



domain. In this context, the operator Aa{w,h), defined in Eq. (72), can be said to be 
the Eulerian counterpart of the operator A-f{w) as is shown by the the following result. 

Theorem 4 (Eulerian and Lagrangian elastic stiffness operators of the immersed do- 



main). With reference to the definitions in Eqs. (72), (78), and (79), we have 



Ac{w,h) ^ Sc^{h)A^{w) and Sa-y{h) ^ M'I^i{h)M:;3 , (80) 

where Sa'y{h)A'y{w) and Ail^i{h)Ai~^ indicate the composition of the operators Say{h) 
and A-y{w) and of the operators M^iih) and , respectively. 



Proof. By the definitions in Eqs. (78) and (79), Vio, h e '3^ and G we have 
,{S^^{h)A^{w),v)^ ^ {A^{w),v{x)\^^ )^^^, 



(81) 



which, using again the definition in Eq. (79), gives 



{Saj{h)Aj{w),v).. 



Pl[w]f^[h]-VM^)L=,^^^,)dV, 



(82) 
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where the second hne of the above equation was obtained by a standard appHcation of 



the chain rule. Comparing the resuh in Eq. (82) with the definition in Eq. (72 1, the first 



of Eqs. (80) foUows. Next, again applying the definition in Eq. (78), yw,h G "3^ and 

(83) 



Vu e %, we have 



which, by the definitions in Eq. (73) and Eq. (75), gives 



J B 



x—s-\-h{s 



dV 



from which we deduce that 



(84) 



(85) 



Since the operator Ad-ys is the Riesz identity between J^/fy and , it is invertible and 
the second of Eqs. (80 1 follows. 



□ 

The operators defined above allow us to formally restate the overall problem described 



by Eqs. (p6|, (54), and (55) as follows: 



Problem 1 (Incompressible fiuid, incompressible solid: dual formulation). Given initial 
conditions Mq G and Wq e iV, for ah t e (0,T) find u{x,t) e p{x,t) £ and 
w{8,t) € ^ such that 

+ 5Mal{'w)u' + 5Afal{w, w' , u)u + 6'Dal{w)u + Sa-y{w)A'y{w) = Fa + Ga{w), 

(86) 

Bp^u = 0, (87) 
M^zw' — M^ii{w)u = 0. 



Remark 10 (Eulerian vs. Lagrangian elastic operators). Referring to Eq. (86), Theorcm|4] 
shows that we could have formulated Problem [l] using the Eulerian elastic operator 
Aa{w^w) instead of the composition Sa-y{w)A^{w). This is because, in the infinite 
dimensional context of our abstract variational formulation, the operators Aa{w,w) 
and Sa-f{w)A^{w) are equivalent. However, as will be shown in Section 3.4 the use of 
Sa-^{w)A^{w) is justified by the fact that this operator lends itself more naturally to the 
derivation of stability estimates that rely solely on the weak form of the velocity coupling 



in Eq. (88). Moreover, anticipating a result discussed in Section |4.4[ it turns out that 



(i) the equivalence between Aa{w,w) and Sa-f{w)A-y{w) fails to hold for the discrete 
version of these operators, and {ii) only the discrete version of Sa^{w)A^{w) can be 
shown to yield a satisfactory semi-discrete stability estimate. 
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3.3. Governing equations: compressible solid 

When the sohd is compressible, incompressibihty must be restricted to the physical 
(as opposed to the extended) fluid domain. In addition, since the stress response in the 
solid is completely determined by the solid's stress constitutive response function, the 
field p contributes to the balance of momentum equation only over the domain D,\ Bt. 
Therefore, for the balance of linear momentum, we write 

/ pf{u — b)-vdv— / pdivvdv+ / Tf -Vx^^dw— / Tg-vda 
Jn Jn Jn JonN 

+ J J{s,t)p{x,t)dWv{x)\^^^^^^^dV 



(89) 



Equation ( 89 1 is identical to Eq. ( 56 ) except for the term appearing as the third line of 



Eq. (89 1. This term can be viewed as a correction to the second term on the first line 
that restricts the contribution of the field p to il\ Bf. 

The restriction of the balance of mass equation to the domain n \ Bf can be written 
as follows: 



q div udv — q div udv = 0. 
n J Bt 



(90) 



To determine the motion of the solid domain, we adopt the same equation presented in 
the case of incompressible solids: 



'w{s, t) — u{x, t) I 



y{s)dV^Q \/y^,ytfY- 



(91) 



Equations ( 89 1-( 91 ) would allow us to determine a unique solution if the field p were 



restricted to the domain VL\Bt. However, our numerical scheme still requires that p be 
defined everywhere in O. To formulate a problem admitting a unique solution for the 
field p € we must sufficiently constraint the behavior of p over Bt- The strategy to 
enforce such a constraint is not unique. In some sense, p can be restricted to fi \ 
by requiring that p — Q over Bt- Another, and perhaps more physically motivated, 
approach is to observe that, for a Newtonian fluid, p represents the mean normal stress 
in the fluid. Therefore, one may choose to constraint the field p in a such a way that it 
represents the mean normal stress everywhere in fi. Since the solid is compressible, its 
mean normal stress is completely determined by the solid's stress constitutive response 
functions. Specifically, letting ps[u,w\ denote the constitutive response function for the 
mean normal stress in the solid, we have 



Ps[U,W\ 



1 

'tTi 



Jl[u]-\ + J-\w]Pl[w]-P[w] 



(92) 
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Therefore, in addition to enforcing Eq. (90), we can enforce the requirement that p — 
Ps[u, w] = over Bt- With this in mind, we replace Eq. (90) with the following equation: 



J qdivudv + J J{s,t)q{x) div u{x,t)\^^^^^^^^dV 
ciJ{s,t)[p{x,t) - C2Ps[u,'w]]q{x)\^ 



dV ^0 Vg e ^, (93) 



where ci > is a constant parameter with dimensions of length times mass over time, 
and where C2 is a dimensionless constant that can take on only the values or 1. For 
C2 = 0, the last term on the left-hand side of Eq. (93) is a (weak) requirement that p — 
over Bf, whereas for C2 = 1, the field p is (weakly) constrained to be equal to the mean 
normal stress in the solid domain and therefore everywhere in Q. 

Remark 11 (True incompressibility vs. near incompressibility) . When a compressible 
solid is immersed in an incompressible fluid, the classes of motions for the solid and 
the fluid, respectively, are not necessarily the same. In this case, problems are typically 
formulated in such a way that the incompatibility is removed by assuming that both the 
fluid and the solid are compressible and then tuning the constitutive parameters of the 



fluid to approximate a nearly incompressible behavior (see, e.g., Blanco et al. 2008a 



Wang et al. 2009). From elasticity it is known that nearly-incompressible models with 



the same incompressible limit behavior may behave differently from one another (see, 
e.g., Levinson and Burgess 1971). For this reasons, the present authors feel that it is 



important to offer a numerical approach to the solution of the problem of a compressible 
solid in an incompressible fluid as its own individual case. 

As was done in the incompressible case, we now reformulate our equations in terms 
of operators defined via duality. Most of the operators defined in the incompressible case 
appear in the formulation of the compressible solid case. Hence, we now define only those 
operators that did not appear in the previous case. Specifically, we define the following 
three operators: 
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,{SBpi{w)u,q)^ 



J[w]q{x) div u{x) 



dV, 



{5Bj,{w)p,v)..^ 



J[w]p{x) diYv{x) 



dV, 



5V, 



I32{ 



^, {S'Pp2iw)p, q)^ 

5£p{u, w,h) e ^ 



J[w]p{x)q{x) 



B 



x—s-\-w{s) 



dV, 



J[h]fl[u] • I + P^H • f[h] q{x)\ 



s+h{6 



(94) 
(95) 
(96) 

dV. 

(97) 



These operators, along with those defined earlier, allow us to formally restate the 
overall problem described by Eqs. (89 1, (93 1, and (91 1 as follows: 
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Problem 2 (Incompressible fluid, compressible solid: dual formulation). Given constant 
coefficients ci > and C2 = V 1, and given initial conditions e 'V and Wq G 'W ^ for 
all t e (0,r) find u{x,t) e 't , p{x,t) G and w{s,t) e ^ such that 

+ 6Mal{'w)u' + 5Afal{w, w' , u)u + 6'Dal{w)u + Sa-^{w)A'y{w) = Fa + Qa{w), 

(98) 

\_Bpi + SBpi (to)] It + ci [5Vi}2{w)p - C28£fs{u, w, w)] = 0, (99) 

Mj3w' - M^i{w)u^O. (100) 



Remark 12 (Complementarity of operators in In Eq. (99), the supports of the terms 
[Bj3i+5Bj3i{w)^ and [SVj32{'w)+C2S£j3{u,w,w)^ are r2\_Bt and St, respectively. That is, 
the supports of the terms in question are complementary subsets of f2. Consequently, the 
terms [Bpi + SBi3i{w)] and [5Vp2i'w) + C2S£p{u,w,w)] are equal to zero individually: 

[Bf3i + SBpi{w)]u ^ and ci[SVf32i'w)p — c2S£i3{u,w,w)] = 0. (101) 

This also implies that the constant ci in Eqs. ( [98| ) and the second of Eqs. ( |101[ ) should 
not be interpreted as a penalization parameter but as a way to ensure that the equations 
are dimensionally correct. 

Problems [l] and [2] can be formally presented in terms of the Hilbert space 3f := 
X ^ X ^V, and 3fo := Yq x ^ x with inner product given by the sum of the inner 
products of the generating spaces. Defining iF 9 ^ := [u,p, wf^ and 3 i/; := [v, q, y]"^, 
then Problems [T] and [2] can be compactly stated as 

Problem 3 (Grouped dual formulation). Given an initial condition ^0 G for all 
t e (0,T) find C(i) e 3f, such that 

(J-(t,e,0,V') =0, VV^GiTo, (102) 

where the full expression of : i-> i^* is defined as in Problem [T] or Problem [2] 

Remark 13 (Initial condition for the pressure). In Problem [sj an initial condition for the 
triple ^0 = [uo,PO}'>^o]'^ is required, just as a matter of compact representation of the 
problem. However, only the initial conditions Uq and Wq are used, since we have no time 
derivative for the pressure which is a Lagrange multiplier for the incompressible part of 
the problem, or completely determined by the solution in the compressible part. 

3.4- Stability of the abstract variational formulation 

The definition of the operators A4ai and A/'qi(m), along with the concept of material 
time derivative for Eulerian fields, yield the following result: 

^,{Maiu',u),^ + .^,{Afai{u)u,u)^^ ^ / ^PfU'^dv. (103) 



n 



Keeping in mind that, for x = s + w{s, t), we have 



d , , \ \ du(x,t) 

U = ■w-U[S + W[S, t),t) — 



dt ' y ' ' 



x—s^'w{s.i) 
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+ V,«(a;,i)| ^ , J'^^^'^K (104) 



and as a straightforward applicatfon of Eq. (1571 in Lemma [S] (with fl replaced by B), 
we have that our definition of the operators (5A^c(i(tt)) and 5j7ai{w , w' , u) is such that 



d 



dt 



B 



yfu^ dv, (105) 



where the above result is due to the fact that we selected w' instead of u in the nonlinear 
advection term of the acceleration of the solid. Therefore, from Eqs. (103) and (1051, 
our formulation is such that 



d 

dt 



^ Pf dv 



n\Bt 

Invoking Eq. (165) in Lemma ([5|, we obtain 

d 



x—s-\-w 



is) 



ipf di> 



n\Bt 



dt 



ipt dv 



fl\Bt 



^Pf u ■ m da. 



dV. (106) 



(107) 



Equations ( 103 ) and ( 105 ) taken together can be written as follows: 



d 

di 



Kdv 



KU ■ m da, 



(108) 



where k — ^pu^. The remaining terms in the stability estimates are the viscous dissipa- 
tive terms 

(109) 



and 



J B 

where the combination of the two yields the dissipative term 



dv. 



{Vaiu.u) + ,{5'Dai{w)u,u)^ / l''[u]-Vudv. 

Jn 



(110) 



(111) 



The final term needed for the derivation of the full energy estimate pertains to the 
time rate of change of the elastic energy. In the proposed immersed method, the coupling 
between the Eulerian and Lagrangian frameworks is embodied by a variety of operators. 
These operators take different forms depending on whether the velocity coupling between 
u and w' is used in its strong form, as in Eq. (39), or in its weak form, as in Eq. (88) 



23 



or ( |100[ ). If we could use directly the strong form of the velocity coupling, namely 
Eg. p9|), we would have that u{s + w{s,t),t) = w'{s,t), and since the chain rule gives 



I x—s-\-w{s) 



F[w], 



(112) 



we would then obtain the usual elastic energy estimates from the definition of the operator 

.^,{Aa.{w,w),u)^^ ^AI^W,e(^f[w])dV. (113) 

However, for solutions u and w' of Problem [T] or [2j Eq. (39) holds only in J^y i that is, 
in its weak form and Ec^. (113) can no longer be obtained as just illustrated. We must 
therefore proceed in a different way. Our starting point is the standard estimate in the 
space for the (fully Lagrangian form of the) stiffness operator Aj{w): 



(114) 



which is valid for w \tv '3^ and w' in J^Y- Using Eq. ([Tm]) and Theorem [1} we can prove 
the following Lemma: 

Lemma 2 (Energy estimate for the immersed elastic operator). Given an elasticity 
operator A^{w), then the operator Sa'y{h)A^{w) satisfies the following energy estimate 
whenever Eq. (881 or Eq. ( 100[ ) are satisfied: 

.^,{Sajiw)A^{w),u)^ - — 



dt 



(115) 



w!if[w])dv. 

Proof. Using Eq. (88 1 or Eq. (100) along with the invcrtibility of A^-y3, the Riesz identity 

(116) 



on we can write 

w' = M:^^M^i{w)u. 
Substituting Eq. ( |116[ ) into Eq. ( |114[ ), we have 

^.{A^{w),M-:^M^iiw)u)^^ - ^ 



W,^{f[w])dV. 



(117) 



Focusing on the left-hand side of Eq. (117), we can writcQ 

^,{Aj{w),m:;^Mji{w)u)^^ = ^,{m^i{w)u,m:;:^Aj{w))^^ 



,{M^,iw)M-^A^iw),u)^, (118) 



where we have applied the definitions in Eqs. (74) and (75) to obtain the last of the 
above expressions. Using the result in Ec^. (118), Eq. (117) can be rewritten as 

^,{M^,{w)M:;^A^iw),u).^ = f Wi{f[w])dV, 



and Eq. (115) follows from the application of Theorem |4] 

^Referring to Eq. \73\ , Ai-,3 is such that ^t(A^^3K;,y) = 
,J^. As a consequence, is such that ^* •^^3^'"*)^^ 
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(119) 
□ 



. (A^73y,if)_-^^ for all w,y e 



for all 



Combining the results in Eqs. (1081, (111), and (115) allows us to state the following 
theorem: 



Theorem 5 (Energy estimate for the abstract variational formulation). Let u, p and w, 
be the solutions of either Problem^ or Problem^ Then the following energy estimate is 
satisfied 



pb ■ udv 



where n = \ pu^ ■ 



u da 



d 
dt 



K dv 



KU ■ m da 



90a 



T'[u]-S/udv+-- W^{F[w])dV, (120) 



Remark 14 (Validity of Theorem [S] in the compressible case). When computing the 
duality between the balance of linear momentum in Eq. (98) and the exact solution u, 



one can use explicitly the properties expressed in Eqs. ( 101 1 to remove the terms involving 



the pressure from the estimate, obtaining formally the same result which is obtained in 
the incompressible case. 

When the external forces and the boundary conditions are identically zero, we obtain 
the following result: 



— / KdV 



and 



dt 



Ku-mda+— I W^dV+ / - S/^udV = 0. (121) 



KdV 



KU ■ m da 



dill, 



dt 



W! dV < 0. 



(122) 



We observe that the presence of the kinetic energy flux (k u ■ m) is due to the fact that 
we have posed our problem over a control volume with mixed boundary conditions. With 
this in mind, Eq. (121) states a classical result: the instantaneous variation of the total 



energy of the system, namely the sum of the kinetic energy, the potential energy, and 
the kinetic energy flux, is equal to the negative of the internal dissipation. Furthermore, 



inequality ( 122 ) implies that for any physically admissible set of constitutive equations for 



the solid and the fluid, the abstract variational formulation we propose is asymptotically 
stable whenever 951 at = 0. 



4. Discrete Formulation 

Spatial Discretization by finite elements 

To approximate the continuous problem, we introduce the decompositions U,h for U, 
and Bh for B into (closed) cells K (triangles or quadrilaterals in 2D, and tetrahedra or 
hexahedra in 3D) such that the usual regularity assumptions are satisfied: 

1. n = \J{K € ^h}, and B = \J{K £ B,,}; 

2. Any two cells K,K' only intersect in common faces, edges, or vertices; 

3. The decomposition U,h matches the decomposition dVt = dflo U S^In- 
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On the decompositions ilh and B^, we consider the finite dimensional subspaces 1^ C 
C ^, and C ^ defined as 



n := e r I e VviK), K€n^}= span{<},^- (123) 

J^h ■■= {ph e ^ I PhiK e 7'Q(if), e O,,} = spanlg^^^fi (124) 
:= e ^ I e Vy{K), KeBh} = span{yj.}^- , (125) 

where Vv{K), Vq{K) and Vy{K) are polynomial spaces of degree ry, rg and ry respec- 
tively on the cells K, and Ny, Nq and iVy are the dimensions of each finite dimensional 
space. 

Our choice of finite dimensional spaces "fh and are included in the pivot spaces 
Jify and respectively, which allow us to use only one discrete space for both u and 
u' and one for w and w' . 

In the examples, we chose the pair and so as to satisfy the inf-sup condition 
for existence, uniqueness, and stability of the approximate solution pertaining to the 



Navier-Stokes component of the problem (see, e.g., Brezzi and Fortin 1991). Note that 



the definitions in Eqs. ( 123 )-( 125 1 imply that the functions in Yh and are continuous 
over flfi and B^, respectively. 

To state the discrete versions of Problems [T] and [2j we first introduce some additional 
notation. Given a discrete functional space, say, fh, one of its elements is identified 
by the column vector of time dependent coefficients uj^{t), j = l,...,iVy, such that 
Uh{x,t) — J2 '^hit)'^hi^) ' where vj^ is the j^^ base element of iji- With a slight abuse 
of notation, we will write M^iUh to mean the multiplication of the column vector it/j by 
the matrix whose elements M^-^ are given by 

where the operator in angle brackets is the one defined earlier. The same is intended for 
all other previously defined operators. 

Remark 15 (Dimension of matrices). The subscript convention that we adopted for the 
continuous operators allows one to determine the dimensions of the matrices and of 
the column vectors involved. For example, the discrete operator Mai is a matrix with 
dimensions Ny x Ny, while the matrix Mji{wfi) has dimensions Ny x Ny. 

Remark 16 (Discrete duality product). With the above notation and due to the linearity 
of the integral operator, we can express duality products in the discrete spaces by simple 
scalar products in M.-^ , where N depends on the dimension of the system at hand. For 
example, given the matrix Mai, then 



{MalUh,Vh)y=Vh- MalUh, (127) 



where the dot-product on the right hand side is the scalar product in . 

For a given choice of ilh and Bh, along with corresponding choices of the finite 
dimensional spaces 1^, ^h, and we reformulate Problem [T] as follows: 
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Problem 4. Given ito £ fh, Wq e for all t € (0,r), find Uh{t) G t^, Ph{t) € ^/i, 
and wi^{t) e such that 

Maiu'^ + Nal{Uh)Uh + DalUh + {B pif Ph 

+ 5Dal{Wh,)ui, + Sa^{Wh)A^{Wh) ^ + Ga{Wh), (128) 

BpiUh - 0, (129) 
M^3w'^ - M^i{wh)uh = 0, (130) 

where u'i^{x,t) = J2[^hi^)]' '^i.^^) '"'/!.(*; ^) — ^'^"^ where the prime 

denotes ordinary differentiation with respect to time. 

Similarly, we reformulate Problem [2] as follows: 

Problem 5. Given constant coefhcients ci > and C2 ~ OVl, and given initial conditions 
Mo e % and Wo € for all t € {0,T) find Uh{x,t) e fh, Ph{x,t) € ^h, and Wh{s,t) e 
^% such that 

MqIM^ + Nal{Uh)Uh + DalUh + [B^i + SBpiiWh)]^ Ph 

+ 5Da,i{wh)uh + 5cj^(u;/j)A^(k;,,) = + Ga(t«/i), (131) 
[S^i + 5Bfii{wh)\ Uh + ci [(5P,32(«^/i)j3h - c2SEi32{uh, Wh, Wh)] = 0, (132) 

M^3w',^~ M^i{wh)uh^O, (133) 

where u'^{x,t) = ■"'h(*'^) — ^i^ii^)]' Vhi^) ^ where the prime 

denotes ordinary differentiation with respect to time. 

In compact notation, Problems |4] and [5] can be cast as semi-discrete problems in the 
space ^ D X -^h ^ 

Problem 6. Given an initial condition G ^h, for all t G {0,T) find ^h{t) G ^h, such 
that 

F{t,^h,eh)^0, (134) 

where 

^^^(t,^^,^) (-^(i,a,e;)>D, i = 0,...,iVv + iVQ+iVy, (135) 



and has the same meaning as in Eq. (102), with ip\ being the basis function for the 
spaces 'Vji, i?/;, or corresponding to the given value of i. 

Theorem 6 (Semi-discrete strong consistency) . The discrete formulations presented in 
Problem^ and Problem^ compactly represented in Problem^ are strongly consistent. 

Proof. The claim follows immediately from observing that, for the exact solution ^ := 
[u,p,w]'^, the equalities 

F\t,^,a-^{Ht,^,e),^'k)=0, t = 0,...,Nv + NQ + NY, (136) 

are satisfied for any conforming approximation, that is, whenever ^ ■^h ^ and 
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4-. 2. Variational velocity coupling 

Earlier in the paper we argued that the use of Dirac-(5 distributions is not a theoretical 
or practical necessity of immersed methods. We now illustrate our implementation of the 
operators embodying the FSI using the standard "infrastructure" of typical FEM codes. 

The operators Mai, Nai{uh), Dai, Bpi, and Fa in Problems |4] and [5] are common 
in variational formulations of the Navier-Stokes problem. We implemented them in a 
standard fashion. The operator M^^ is the mass matrix of the space J% and, again, its 
implementation is standard. The non-standard operators in our formulation are those 
with a nonlinear parametric dependence on the field w (the motion of the solid). We now 
discuss the construction of the matrix M^i{w) (corresponding to the operator defined 



in Eq. (74)) which is responsible for a successful couphng of velocities between the fluid 
and the solid domain. 

For convenience, we recall that the entries of the matrix Mryi{wh) are given by: 

= <^sj <(^)L..+.,(,,)-yU«)dK (137) 



The construction of M^i{wh) requires that we compute the integral in Eq. (137). As 
is typical in FEM, this is done by summing the contributions due to each cell K of 
the triangulation Bh. These contributions are computed using quadrature rules with 
Nq points. That is, the contribution of an individual cell is computed by summing the 
value of the products of the integrand at the quadrature points times the corresponding 
quadrature weight. The integrand consists of the functions y\{s), whose support is 
defined over the triangulation of Bh, and of the functions v\{x) (with x = s + Wh{s,t)) 
whose support is instead defined over the triangulation flh- Operationally, we perform 
this calculation as follows. First we determine the position of the quadrature points of 
the solid element, both relative to the reference unit clement and relative to the global 
coordinate system adopted for the calculation, through the mappings: 

SK:K:^[0,lf^KeBh, (138) 
I + Wh-.K ^ solid ceU. (139) 

Next, the global coordinates of the quadrature points (obtained through the mappings 
in Eqs. (138) and ( 139[ )) are passed to a search algorithm that identifies the fluid cells 



in Qh containing the points in question, at which the functions vj^ are evaluated. The 
outcome of this operation is sketched in Fig. [2] where, as a way of example, we show the 
image (under the motion s + Wh{s, t)) of a cell of Bh straddling four cells of denoted 
fluid cells A-D. The quadrature points over the solid cell are denoted by fllled circles. 
The contribution to the integral in Eq. ( |137 ) due to the solid cell is then computed by 



summing the partial contributions corresponding to each of the fluid cells intersecting 
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• : quadrature point of the solid cell 



fluid cell C 




fluid cell D 
/ — solid cell 


i 

fluid cell A 




u 

fluid cell B 



Figure 2: Cells denote as A-D represent a four-cell patch of the triangulation of the fluid domain. The 
cell denoted as "solid cell" represents a cell of the triangulation of the immersed solid domain that is 
contained in the union of cells A-D of the fluid domain. The filled dots represent the quadrature points 
of the quadrature rule adopted to carry out integration over the cells of the immersed domain. 



the solid cell in question: 



NK.g 



yUsK,qpK.q, (140) 



where we denoted with sx < 



the transformation of the g-th quadrature point under the 
the corresponding quadrature weight. 



mapping sk, defined in Eq. (138), and with ijJK,q 
In general, the number of quadrature points corresponding to each partial contribution 
varies. The implementation of an efficient search algorithm responsible for identifying 
the fluid cells that deflne the partition of an individual solid cell is the only technically 
challenging part of the proposed immersed method. However, several standard techniques 
are available to deal with this task (see, e.g 



Thompson et al. 


1998 


de Berg et al. 


2008 



Once the fluid cells containing the quadrature points of a given solid cell are found, we 
determine the value of vj^ at the quadrature points using the interpolation infrastructure 
inherent in the flnite element representation of fields defined over flh- Our finite element 



code was developed using the finite element library deal. II (Bangerth et al., 20071 



which provides built-in facilities to carry out precisely this type of calculation. 



4-.3. Variational force coupling 

The coupling between the Eulerian and Lagrangian frameworks is embodied by the 
operator Sa-yiw). The discrete version of this operator is constructed using the discrete 
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versions of the operators A4ji{w) and Ai^s, and Theorem |4j 



(141) 



where M^^ is the usual mass matrix for the space and M^i(wh) is the transpose of 



the coupling matrix discussed in Section [4~2 
As we observed in Remark 



10 



the operators Aa{h,w) and Sa-y{h)Aj{w) are equiv- 
alent in the abstract variational formulation. At the discrete level, however, this is no 
longer the case. When approximating Aj{w), one needs to integrate terms that contain 
the gradient of the basis functions yj^ of the space "S^h- By contrast, the approximation of 
Aa{h, w) requires the evaluation of the gradients of the basis functions w)^, in the space 
fh, under the map s + Wh{s,t). In general, we have that 

Aa{hh,Why := ^.{Aa{hh,Wh),v{)_^ ^ [M^^{hh)M-^ A^(wh)y =: {So,^{hh)A^(wh)y 

(142) 

and no equivalence can be shown in the discrete space between the two discrete operators 



in Eq. (142). In principle one could use in the discretization the natural definition of 



the discrete operator Aa{hh, Wh). However, only the right hand side of Eq. ( 142 ) can be 
shown to satisfy a discrete energy estimate equal to that in Eq. (113): 



Theorem 7 (Discrete energy estimate for immersed elastic operator). Given an elasticity 
operator A-yiw) and its discrete counterpart Ary(wh), then the discrete operator 



S^^{hh)A^{wh) := M^,{hh)M^^A^{wh) 



(143) 



satisfies the following semi-discrete energy estimate, whenever Eq. ( 130 ) or Eq. { 133 1 are 
satisfied: 

{Sa,^{wi^,)A^{wi^,)) .Uh^^^J^ Wiif[wi^,]) dV. (144) 

Proof. The proof follows closely that of Lemma [2] If we take the scalar product of the 
semi-discrete version of the velocity coupling Eq. ( |l30[ ) with the term (^M~^Aj{wfi)y we 
obtain 

{M-^^A^{wh)) ■ M^3w'^ - {M:;^^A^{wh)) ■ M^i{wh)uh 

= A^iwh) ■ w'^ - {M^^iwh)M-:,'A^{wh)) ■ = 0. (145) 

The discrete estimate deriving from Eq. ( |114[ ) then gives immediately the semi-discrete 
estimate of Eq. (|144| . □ 



Remark 17 (Spread operator). Using the approximation strategy described in Eq. (142), 
the only operator that couples directly the Eulerian and the Lagrangian framework is 
A4.yi{w), whose implementation details have been discussed in Section 4.2 Notice that it 
is essential to use the same operator in the momentum conservation equation (specifically, 
its adjoint) to obtain the discrete stability estimate in Eq. (144 1. In the IBM literature, 
the adjoint of M-^i^w) is also known as the spread operator, because of its role in 
distributing the forces due to the elastic deformation of the immersed domain to the 
underlying fluid domain. 
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4.4- Semi discrete stability estimates 

Repeating all passages from Eq. (103) to Eq. (1201 in the discrete space iji, we 
wish we could show semi-discrete stability estimates equivalent to those of the abstract 
variational formulation. Unfortunately, contrary to what can be done in the continuous 
case, in the discrete problem we cannot invoke Eq. ( 165 ) in Lemma ^ to say 



n\Bt 



\pi ul Av = 



d 
dt 



\pf u\dv 



n\Bt 



5/9f ul Uh 



m da. 



(146) 



This is a well known fact in the discretization of the Navier-Stokes equations. That is, 
there are stability issues related with the non-linear transport term Mai(u) defined in 
Eq. (65 1. These stability issues originate from the fact that the approximation of Eq. ([T]) 
in the numerical scheme is not satisfied pointwise. And it is this fact that prevents us 
from a direct application of Theorems [2] and [3j With this in mind, we also know that 
stabilization techniques for the operator in question exist that lead to stable formulations 

1982j). Therefore, in the present 



N 



(see, e.g., Heywood and Rannacher 



when dVL 

paper we will limit ourselves to appealing to such stabilization techniques and assume 



that Eq. ( 146 ) is satisfied also at the discrete level. 



Theorem 8 (Semi discrete energy estimate). Let Uh, Ph o,nd Wh, be the discrete solutions 
of either Problem^ or Problem^ Assuming that a stabilized non-linear term Nai{uh) 
is used, such that Eq. ( |146[ ) is satisfied, then the following semi discrete energy estimate 
is satisfied 



pb ■ Uhdv 



where = \pu\. 



Tg ■ Uhda^ 



an I, 



d_ 
d< 



Kh dv 



KUfi ■ m da 



80.1. 



^T"K].VM^dz;+- / WI{V[wn])dV, (147) 



4.5. Time discretization 



Equation ( 134 ) represents a system of nonlinear differential algebraic equations (DAE), 



which we solve using the IDA package of the SUNDIALS OpenSource library ( Hindmarsh 



et al. 


2005 


I. 


et al. 


200511^ 



Hindmarsh 



The integration method in IDA is variable-order, variable-coefficient BDF 
[backward difference formula], in fixed- leading-coefficient form. The method 
order ranges from 1 to 5, with the BDF of order q given by the multistep 
formula 



9 

E 

i=0 



(148) 



^We quoted directly from the SUNDIALS documentation. However, we adjusted the notation so as 
to be consistent with ours and we numbered equations according to their order in this paper. 
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where ^„ and are the computed approximations to ^(i„) and re- 
spectively, and the step size is /i„ = t„ — The coefficients a„^i are 
uniquely determined by the order q, and the history of the step sizes. The 
apphcation of the BDF [in Eq.] (|148[) to the DAE system [in Eq.] (|134[) results 



in a nonhnear algebraic system to be solved at each step: 



(149) 



Regardless of the method options, the solution of the nonlinear system [in 
Eq.] ( 149 1 is accomplished with some form of Newton iteration. This leads 



to a linear system for each Newton correction, of the form 

n,m+ 1 ^ S,n,m] — ^G{^n,m), 



(150) 



where S,n,m is the mth approximation to S,m- Here J is some approximation 
to the system Jacobian 



J = 



dG 



dF 

9e 



dF 



(151) 



where a = ctn.o/hn- The scalar a changes whenever the step size or method 
order changes. 

In our finite element implementation, we assemble the residual G{^n,m) at each Newton 



correction, and let the Sacado package of the Trilinos library (Bartlett et al. 2006 Gay 



1991 Heroux et al. 2003) compute the Jacobian in Eq. (151). The detailed procedure 



used in our code to compute the Jacobian through Sacado was taken almost verbatim 



from the tutorial program step-33 of the deal. II library (Bangerth et al. [1998-2006 1 



The final system is solved using a preconditioned GMRES iterative method (see, e.g. 



Golub and Van Loan 1996). 



5. Numerics 



We present a numerical experiment designed to test the main characteristics of the 
proposed immersed method and those elements that distinguish it from other methods 
in the literature. Specifically, we consider the case of a solid with mass density different 
from that of the fluid. The solid is assumed to be compressible and viscoelastic. The 
dynamic viscosity of the solid is taken to be twice that of the fluid and the elastic part 
of the behavior was chosen to be a compressible neo-Hookean material. The fluid is 
modeled as truly incompressible, as opposed to nearly incompressible. 

5.1. Discretization 

The approximation spaces we used in our simulations are the piecewise bi-quadratic 
spaces of continuous vector functions over VL and over B for the approximations of the 
velocity field Uh and of the displacement field Wh (usually referred to as the continuous 

space), and the piecewise discontinuous linear space over i7 for the approximation 
of the pressure field p. 
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The — pair of spaces is known to satisfy the inf-sup condition for the approx- 



imation of the Navier-Stokes part of our equations (see, e.g., Brezzi and Fortin, 19911, 
while the choice of the space for the displacement variable Wt is a natural choice, 
given the underlying velocity field Uh- With this choice of spaces, Eqs. (133) and (130) 



can be satisfied exactly when the solid and the fiuid meshes are matching. 

One of the advantages of immersed methods is the possibility to select the meshes 
over the fluid and the solid domains independently. However, accuracy issues may arise 
if the mesh over the solid domain is not sufficiently refined relative to that for the fluid 
domain. It has been observed (see, e.g., Peskin 2002) that a reasonable choice is to take 
the mapped solid mesh size hs to be at least one half of the fluid mesh size /if. This 
choice finds its justification in the approximation properties of both the velocity and the 
force coupling schemes presented in Section |4.2| and Section |4.3| It is essential for the 



success of immersed methods that the integrals presented in Eq. ( 140 ) be approximated 
as accurately as possible. Independently on the choice of approximating spaces, there 
will be errors in the approximation of these integrals due to the non-matching nature of 
the fiuid and solid meshes. If one uses a fixed number of quadrature points (as in our 
case), reducing hs while maintaing hi constant increases the accuracy of those integrals 
up to the point in which one element of the solid mesh is entirely contained in an element 
of the fluid mesh. Further reduction of the solid mesh size beyond this point is not useful, 
since it only increases the computational cost, without adding accuracy to the method, 
which is bounded anyway by the fluid mesh size /if. 

The choice /ig ~ | /if is a reasonable compromise, for which most of the solid elements 
are fully contained in a fluid element, and each solid element spans at most four elements 
of the fluid mesh. At run time, whenever a solid mesh element is distorted to span more 
than four fluid mesh elements, the element in question should be refined to increase the 
accuracy of the method. Currently, such tests are not implemented in our code, and we 
select a slightly finer solid mesh to prevent distortion from causing a drift in the accuracy 
of the method. 

An alternative solution is to use adaptive quadrature rules in the approximation of the 
integrals in Eq. (140), as done, for example, in Griffith and Luo (2012). This approach 
allows one to choose hs independently from /if, and it works effectively even in the case 
where the solid cell spans several fluid cells. Conservation of mass in this case may be an 
issue, since the details at which the fluid evolves in the background may not be captured 



accurately enough by the solid mesh (see Griffith ( 2012 ) for detailed discussion on volume 



conservation in Immersed Boundary Methods). 

5.2. Constitutive settings 

We present a simple numerical example concerning a two-dimensional rubber disk, 
modeled by a viscoelastic compressible material, where the elastic part of the behavior is 
that of a compressible neo-Hookean material. The disk is pre-deformed with a uniform 
compression that changes its diameter to a fraction of its diameter in the reference 
configuration, and then it is released from rest in a two-dimensional box containing 
a fluid, also at rest. The dynamic viscosity and mass density of the fluid are /if — 
10~^ kg/(m^ -s) and density pf = Ikg/m^, respectively. On the top side of the box 
we impose homogeneous Neumann boundary conditions, to allow the fluid to enter and 
exit the box, while on the other three sides we impose homogeneous Dirichlet boundary 
conditions. 
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The reference configuration of the solid is a disk of diameter = 0.125 m, centered 
at the origin. Its initial displacement field is given by 



Wo 



-0.3si 
-0.3s2 



0.6 m 
0.4 m 



(152) 



where si S2, expressed in meters, denote the coordinates of points in B relative to the 
chosen Cartesian coordinate system 
function of the solid is T, = T; 



Referring to Eq. (17), the constitutive response 
tr with 



G 



j-2u/il-2u)^-T] Jl=2^,,D, 



(153) 



where G = 20Pa-m, v = 0.3, /Xg = 2 x 10~^ kg/(m^ -s). The mass density of the solid in 
the reference configuration is — O.Spf . We add a constant external body force density 
(gravity) directed downwards: 





-lOmVs^ 



(154) 



The initial deformation of the disk is such that its density (in the deformed state) exceeds 
that of the surrounding fluid. Under these conditions, the disk would sink. However, as 
soon as the disk is released, the disk will expand rapidly and regain a size such that the 
disk will start floating almost from the start of the motion. 

5.3. Results 

In Figs. |3] and |4] we show the evolution of the pressure and of the velocity fields, 
respectively. The plots in Fig. [3] show the mean normal stress in both the fluid and the 
solid, as defined in Eq. (92 1. 



In the figures, three different phases can be recognized: 

1. expansion, for < t < 0.3 s; 

2. contraction and ascension, for 0.3s < t < 2s; 

3. expansion and rising, for 2 s < t < 3 s. 

In phase one, the disk tries to reach an equilibrium state of deformation by quickly 



expanding and pushing the surrounding fluid. Fig. 4(a) shows a snapshot of the velocity 



fleld in this phase: an outflow is present at the top of the box, and the radial velocity in 
the solid shows the role of compressibility in the solid constitutive behavior. Quantitative 



measurements can be inferred from Fig. 5(a) displaying the plot of the total instantaneous 



flux through the Neumann part of the boundary: 



F := j u ■ n ds. 



(155) 



In phase two, the solid has reached a positive buoyancy status, and starts moving 
towards the top of the box. In this phase, the vertical position of the center of mass of the 
disk (Fig. |5(b)] ), increases in an approximately quadratic manner with time. The area of 
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(a) t = (b) t = .5 s 




(e) t = 2s (f) * = 2.5s 



Figure 3: Pressure evolution. 
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(e) t = 2s (f) * = 2-5s 

Figure 4: Velocity evolution. 
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the disk bounces back due to inertial effects (lighter hue in Fig. 5(c) ), and in phase three, 
it grows again, both for a bouncing effect and because of the reduced pressure applied 
on the surface of the disk itself. 

By tracking the vertical location of the center of mass of the disk (Fig. 5(b)), we 
observe that the dynamics of the expansion phase are rather fast, and the disk expands to 
a positive buoyancy state while remaining substantially still. We monitor the consistency 
of the method by computing the integral in time of the total flux of fluid through the 
top side, which should equate the area change of the disk (Fig. 5(c) I: 

SA[:= [ F(T)dT= [ [ u{T)-ndsdTK [ J[w{t)]dv- [ J[wQ]dv := SA^. (156) 
Jo Jo Jfn Jb Jb 

While the consistency of the method in phase one is quite good, the fluid flux does not 
seem to compensate accurately for the small changes in area due to the bouncing of the 
disk in the remaining two phases. This lack of accuracy is likely due to the combination 
of the errors in the approximation of the divergence free constraint in the fluid, and to 
the errors in the computation of the velocity coupling between the fluid and the solid. 



6. Summary and Conclusions 



We presented a fully variational formulation of an immersed method for the solution 
of FSI problems. Like other immersed methods, ours is based on the idea of keeping 
independent discretizations for the fluid and for the solid domains. The fluid is treated 
in its natural Eulerian framework, while the solid is modeled using a Lagrangian strategy. 



Most of the implementations of immersed methods refer to the pioneering work oflPeskin 



( 19771, in which a clever reformulation of the continuous coupling between the fluid and 



the solid domain allows one to construct projection operators between the Lagrangian 
and the Eulerian framework based on approximated Dirac-(5 distributions. While the 
necessity to introduce approximated Divac-d distributions is strongly connected to the 
particular approximation strategy chosen to discretize the continuous problem (FD in 
the IBM), its use has propagated also in the Finite Element community (see, e.g., Zhang 
et al.[ [2004 ). A variational approach that removed the necessity to approximate the 



Dirac-(5 distribution has been proposed in 
inlHeltail p006|; IBofR et al.l (pOOTl |2008|; 



Heltai (2008) 



BofR and Gastaldi (2003), and later extended 



The formulation we presented extends that of Boffi et al. ( 2008 ) to general elasticity 
problems in which the solid and the fluid can have different mass densities. In addition, 
the constitutive response function for the solid can be either compressible or incompress- 
ible, and either viscoelastic of differential type or purely elastic. The abstract variational 
formulation we proposed is shown to yield energy estimates that are formally identical 
to those in the classical context of continuum mechanics. The numerical approximation 
we proposed is strongly consistent and stable, with semi-discrete energy estimate that 
are formally identical to those of the abstract variational formulation and therefore to 
those in the classical context of continuum mechanics. 

We discussed in details the algorithmic strategies for an efficient implementation 
of the proposed method, and we showed how standard implementations of the finite 
element method along with some appropriate search algorithm for the determination of 
the element containing a given point are enough to implement the proposed formulation. 

37 





(c) Area exchange compaj-ison. (d) Area exchange error. 

Figure 5: Instantaneous flux, position, and area exchange. 
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A simple numerical experiment was used to test the novel characteristics of our 
method. While the results are promising, some work is still necessary to ensure bet- 
ter conservation properties of the method. 



Acknowledgements 

The research leading to these results has received specific funding under the "Young 
SISSA Scientists' Research Projects" scheme 2011-2012, promoted by the International 
School for Advanced Studies (SISSA), Trieste, Italy. 

The final version of this paper was much improved thanks to suggestions made by 
the anonymous reviewers. 



Appendix: Proof of Theorem [2] 



The proof of Theorem [2] is presented at the end of this appendix and is preceded by 
some useful intermediate results. 

The application of Theorem [l] to the domains ft and Bt defined in Section |2.1| yields 
the following results: 

Lemma 3 (Transport theorem for fixed control volumes and for physical bodies). Let 

and Bt, with outward unit normals m and n, respectively, be the domains defined in 
Section 2.1 Let <j){x, t) and ^(a;, t) he a smooth field defined over fl and Bt, respectively. 



Then we have 



and 



dt 



[ c^ix,t)dv^ [ 
Jn Jn 



d(j){x, t) 
dt 



dv 



^ f ^{x,t)dv^ f ^^^^/^ dv+ f ^{x,t)u-nda. 



(157) 
(158) 



Proof of Lemma^ The results in Eqs. (1571 and (1581 are well known. The proof is 
presented simply to facilitate the discussion of subsequent results. Since is a fixed 



control volume, its boundary is time independent. Hence, Eq. ( 157) follows directly from 
Eq. (22) when we let u — 0. Next, we observe that Bt is a time-dependent domain such 



that the velocity field on the boundary of Bt coincides with the material velocity field u. 
Hence, Eq. ( |158 ) follows from Eq. (22) when we set u = u. □ 



Lemma 4 (Transport theorem for n\Bt). Let and Bt, with outward unit normals m 

Let 4>{x,t) he a smooth field 



and n, respectively, he the domains defined in Section 2.1 
defined over the domain D,\ Bt. Th- 



en we nave 



d_ 
dt 



(a;, t) dv 



n\Bt 



n\Bt 



d(j){x, t) 
dt 



dv 



[x, t)u ■ n da. 



(159) 



dBt 
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Proof of Lemma We observe that 



d{n\Bt) =dn\JdBt 



(160) 



The unit normals outward relative to n\Bt on dil and dBt are m and — n, respectively. 
FinaUy, the velocity field ol dil is null whereas the velocity field on dBt is equal to the 
material velocity field u (on dBt). Then the result in Eq. (159) follows directly from 
Eq. (22 1 by setting i/ = on d^l and = tt on dBt- 



□ 



The coordinated application of the transport theorems with the balance of mass and 
the concept of material time derivative yields results that are useful in the derivation 
of energy estimates. If <j>{x,t) is the Eulerian description of a scalar-valued physical 
quantity, then the material time derivative of (f> is 



dt 



+ grad (f) ■ u, 



(161) 



where we recall that u{x,t) is the Eulerian description of the material velocity field. We 
now consider the case in which (f> is a density per unit volume with corresponding density 
per unit mass V'j so that (j){x,t) = p{x,t)ip{x,t), where p{x,t) is Eulerian description of 
the mass density distribution. Then, Eq. (161) gives 



pip + pijj 



djptp) 
dt 



+ gia.d{pip) ■ u 



dt 



pip + pip — grad(pV) • (162) 



Using Eq. ([T]) and recalling that div(0tt) = grad0 • u + (pdivu, the last of Eqs. (162) 
becomes 

^^^piP-dW{piPu). (163) 

Lemma 5 (Transport theorems for densities per unit mass). Let f2 and Bt be the domains 
defined in Section 2.1 Let ipBt{x,t) and ipn\Bti^T^) Eulerian descriptions of 

sufficiently smooth scalar-valued physical density per unit mass defined over Bt, and 
^\ Bt, respectively. Then, Theorem^ and the principle of balance of mass in Eq. ^ 
imply 



d 
dt 



piJjBt dv ^ pipBt dv 



Bt 



d 

dt 



pA2\Bt<^v+ pip^\BtU-mda^ pipn\Bt'^v. 
n\Bt Jon Jn\Bt 



(164) 
(165) 



Proof of Lemma^ Equation (164) is a well known result that can be found in many 
textbooks (see, e.g., Gurtin et al. 2010). It is obtained by setting ^ — pipBt in Eq. (158) 



and then using Eq. ( 163 ) along with the divergence theorem. This same strategy can be 



used to obtain Eq. ( 165 ), that is, substituting pipn\Bt in place of ( 
Eq. ( 163 ) along with the divergence theorem. 



in Eq. ( 159 ) and then 
□ 



Proof of Theorem^ Equation (23) is obtained by summing Eqs. (164) and (165), where 
ipBt and ipa\Bt ^^'^ taken as the restrictions of ip to Bt and fl\ Bt, respectively. □ 
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